Example of Calibration¶

Import stuff

# Default

budget = 8192
num_workers = 256
maximum_jobs = 64

import warnings
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

from dask.distributed import Client
from dask_jobqueue import SLURMCluster as Cluster

import nevergrad as ng
from mhpc_project.matsch_b2 import CalibrationModel, Variables, Loss, NGO
from geotopy.utils import date_parser, comparison_plot, ProgressBar, ParametersLogger
from geotopy.measures import KGE

Load observations and apply some filters

observations = pd.read_csv('../data/Matsch B2/observations/obs.csv',
                           na_values=['-9999', '-99.99'],
                           usecols=[0, 2],
                           parse_dates=[0],
                           date_parser=date_parser,
                           index_col=0,
                           squeeze=True)
observations.index.rename('datetime', inplace=True)

selection = observations.index.map(lambda x: not (x < date_parser('01/03/2011 23:00') or x.month in [12,1,2]))

observations = observations[selection]

Load observations, model, loss function and calibration strategy

model = CalibrationModel('../data/Matsch B2/geotop', run_args={'timeout': 180})
variables = Variables('../data/Matsch B2/variables/testbed.csv')
measure = KGE(observations)
loss = Loss(model, variables, measure)
calibration = NGO(loss, budget=budget, num_workers=num_workers)

Evaluate the model and plot the risults

simulation = model()
print(f"Before optimization loss is {measure(simulation)}")
comparison_plot(observations,
                simulation,
                desc='Soil moisture content @ 5cm',
                rel=True)
plt.show()

Before optimization loss is 0.7071019548837151

Launch Dask cluster using dask-jobqueue

cluster = Cluster()
cluster.adapt(maximum_jobs=maximum_jobs)

client = Client(cluster)

Run calibration

logger = ParametersLogger(loss.massage)

calibration.optimizer.register_callback('tell', logger)

calibration.optimizer.register_callback('tell', ProgressBar(total=calibration.optimizer.budget))

with warnings.catch_warnings():
    
    warnings.simplefilter("ignore")
    recommendation = calibration(executor=client)

(128_w,256)-aCMA-ES (mu_w=66.9,w_1=3%) in dimension 10 (seed=<module 'time' (built-in)>, Fri Nov 27 19:23:25 2020)

Plot results and show table of calibrated values

simulation = model(**recommendation)
print(f"After optimization loss is {measure(simulation)}")
comparison_plot(observations,
                simulation,
                desc='Soil moisture content @ 5cm',
                rel=True)
plt.show()

After optimization loss is 0.2401898050231394

calibration.to_dataframe(recommendation)

	lower	upper	recommendation
NormalHydrConductivity	0.00001	1.0	0.013477
ThetaRes	0.00000	1.0	0.189842
WiltingPoint	0.00000	1.0	0.013288
FieldCapacity	0.00000	1.0	0.607932
ThetaSat	0.00000	1.0	0.647526
NVanGenuchten	1.01000	2.5	1.147011
AlphaVanGenuchten	0.00001	1.0	0.000860
LSAI	0.00000	5.0	1.042750
RootDepth	0.00000	1500.0	368.257241
SoilRoughness	0.10000	100.0	0.718859

logger.parallel_coordinate_plot()

HiPlot

Loading HiPlot...

from bokeh.plotting import figure, show
from bokeh.io import output_notebook

output_notebook()

Loading BokehJS ...

p = figure(x_axis_type='datetime', tools='xpan, xwheel_zoom', active_scroll='xwheel_zoom')
p.line(observations.index.to_numpy(), observations.to_numpy(), line_color='blue', line_width=2, legend_label='observations')
p.line(simulation.index.to_numpy(), simulation.to_numpy(), line_color='red', line_width=2, legend_label='simulation')
show(p)

An empirical study on the calibration of the GEOtop model using supercomputers

Example of Calibration¶